import numpy as np
from MyDeFunc import MyDeFunc
import matplotlib.pyplot as plt
import os
from datetime import datetime


# 定义常数
a_coal=0.85 
a_grid_e=a_coal #定义电网NCI
a_chp_ec=0.50
a_gas=1.96
a_gb=0.45
c_h=[0.094,0.030,0.016,0.015,0.025,0.039] #设备h的运行费用
c_gas=1.964 #天然气的碳排放
P_egrid_max=5000 #电网买，卖电最大功率
P_ggrid_max=5000 #气网买，卖气最大功率
LHV_ng=35.6 #天然气低热值

P_chp_max=1500 #CHP容量，功率限制都是0-最大功率 单位kW
P_chp_ef=0.92
P_chp_cost=0.094 #运行消耗，单位CNY/kW
P_gb_max=1000
P_gb_ef=0.43
P_gb_cost=0.030
P_eb_max=2000 
P_eb_ef=0.98
P_db_cost=0.016
P_p2g_max=1000 
P_p2g_ef=0.64
P_p2g_cost=0.015
P_wt_cap=1500

P_ees_cap=1000 # EES存储容量，单位kW
P_ees_cap_t=np.full(24,P_ees_cap)
SOC_EES_T=[]
n_ees_ch=0.90
n_ees_dis=0.90
P_ees_max=1000 #EES充放电功率
P_tes_cap=500
P_tes_cap_t=np.full(24,P_tes_cap)
SOC_TES_T=[]
n_tes_ch=0.95
n_tes_dis=0.90
P_tes_max=500
P_ges_cap=500 #单位立方米
P_ges_cap_t=np.full(24,P_ges_cap)
SOC_GES_T=[]
n_ges_ch=0.95
n_ges_dis=0.95
P_ges_max=300

# 设备转换效率
n_t_t2e=1.47 # CHP的热电联产比
n_t_e2t=0.6 # EB的电热转换比率
n_p2g=0.6 # P2G的电气转换效率
n_gb=0.6
n_chp=0.6

#不同设备功率
P_eload_t=np.array([3800,2400,1050,2900,2850,3800,2100,2300,3300,2300,2100,2300,
                    3000,2400,2700,2600,2500,2100,2200,1600,1700,2300,2200,1600]) #用户电负荷
P_tload_t=np.array([1200,1300,1200,1250,2150,2550,2700,3100,1800,1600,1950,1700,
                    1500,1300,1350,1450,2400,2300,1800,1950,2100,1800,1300,1200]) #用户热负荷
P_gload_t=np.array([750,720,700,720,720,850,1100,1300,1400,1150,1130,1120,
                    1050,1170,1050,1080,1500,1300,1350,1300,1100,1000,950,900]) #用户气负荷
P_pv1_t=np.array([0,0,0,0,0,0,0,50,120,190,210,210,205,200,200,190,170,140,20,0,0,0,0,0])
P_pv2_t=np.array([0,0,0,0,0,0,0,0 ,0  ,100,150,200,220,220,200,200,180,150,50,0,0,0,0,0])
P_pv3_t=np.array([0,0,0,0,0,0,0,0 ,50 ,110,120,130,140,170,160,150,120,10 ,0 ,0,0,0,0,0])
P_pv4_t=np.array([0,0,0,0,0,0,0,0 ,0  ,20 ,50 ,60 ,90 ,100,100,100,110,115,50,0,0,0,0,0])
P_pv5_t=np.array([0,0,0,0,0,0,0,50,120,210,270,300,320,330,310,280,200,110,30,0,0,0,0,0])
P_pv_t=P_pv1_t+P_pv2_t+P_pv3_t+P_pv4_t+P_pv5_t
# P_eload_t_net=P_eload_t-P_pv_t #电网功率情况
# P_chp_t=[] #CHP的电功率，热功率为n_t_t2e*P_chp_t
P_wt_t=np.array([0,0,0,0,0,0,0,500,1200,1900,2100,2100,2050,2000,2000,1900,1700,1400,200,0,0,0,0,0])

k=0.04 #碳税CNY/kg
c_esp_t_e=[1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.75,0.86,0.49,0.45,0.43,0.45,0.42,0.50,
           0.59,0.79,0.99,0.99,0.99,0.99,0.99,0.99] # 不同时刻EPS定义的初始电价
c_esp_t_ee=[] #不同时刻产消者出售给ESP电价
c_esp_t_t=[0.46,0.46,0.46,0.46,0.46,0.47,0.49,0.50,0.50,0.47,0.49,0.46,0.46,0.46,0.46,0.46,
           0.46,0.47,0.50,0.50,0.50,0.46,0.46,0.47] # 不同时刻ESP定义的初始热价
c_esp_t_g=[3.45,3.50,3.50,3.50,3.45,3.50,3.50,3.50,3.45,3.45,3.45,3.50,3.48,3.50,3.48,3.48,
           3.50,3.50,3.50,3.50,3.50,3.50,3.45,3.45] # 不同时刻ESP定义的初始气价
c_ebuy_t=np.full(24,0.4) #不同时刻电网的购电价格
c_esell_t=np.full(24,1.0) #不同时刻电网的售电价格
c_gsell_t=np.full(24,0.5) #不同时刻气网的售气价格
c_tsell_t=np.full(24,3.5) #不同时刻热网的售热价格


def verifyEES(x):
  target=np.zeros(24)
  for i in range(1,24):
    target[i]=x[13][i]-x[13][i-1]-(n_ees_ch*x[7][i]-x[8][i]/n_ees_dis)
  return target

def verifyTES(x):
  target=np.zeros(24)
  for i in range(1,24):
    target[i]=x[14][i]-x[14][i-1]-(n_tes_ch*x[9][i]-x[10][i]/n_tes_dis)
  return target

def verifyGES(x):
  target=np.zeros(24)
  for i in range(1,24):
    target[i]=x[15][i]-x[15][i-1]-(n_ges_ch*x[11][i]-x[12][i]/n_ges_dis)
  return target

#约束
constraint_eq=[
  lambda x:P_eload_t-x[0]-P_wt_t-x[6]-x[8]+x[7]+x[5]+x[2]+x[3], #31
  lambda x:P_tload_t-n_t_t2e*x[0]-n_t_e2t*x[2]-x[4]-x[10]+x[9], #32
  lambda x:P_gload_t-x[1]-n_p2g*x[3]-x[12]+x[11]+x[0]/(n_chp*LHV_ng)+x[4]/(n_gb*LHV_ng), #33
  lambda x:verifyEES(x), 
  lambda x:verifyTES(x), #36
  lambda x:verifyGES(x),                  # 电池等式约束
  lambda x:x[5]*x[6],
  lambda x:x[7]*x[8],
  lambda x:x[9]*x[10], #35
  lambda x:x[11]*x[12],
]
L=np.full(24,1)

constraint_ueq=[
  lambda x:-x[6],
  lambda x:-x[5],
  lambda x:-x[1],
  lambda x:-x[7],
  lambda x:-x[8],
  lambda x:-x[9],
  lambda x:-x[10],
  lambda x:-x[11],
  lambda x:-x[12],
  lambda x:x[6]-P_egrid_max,
  lambda x:x[5]-P_egrid_max,
  lambda x:x[1]-P_ggrid_max,
  lambda x:x[7]-P_ees_max,
  lambda x:x[8]-P_ees_max,
  lambda x:x[9]-P_tes_max,
  lambda x:x[10]-P_tes_max,
  lambda x:x[11]-P_ges_max,
  lambda x:x[12]-P_ges_max,
  # lambda x:x[2]+x[3]-x[0]-P_wt_t-x[10],
  # lambda x:x[8]/n_ees_dis-x[13], #放电功率小于当前电量
  # lambda x:x[7]/n_ees_ch-P_ees_cap_t+x[13], #充电功率小于当前空闲
  # lambda x:x[10]/n_tes_dis-x[14], #放电功率小于当前电量
  # lambda x:x[9]/n_tes_ch-P_tes_cap_t+x[14], #充电功率小于当前空闲
  # lambda x:x[12]/n_ges_dis-x[15], #放电功率小于当前电量
  # lambda x:x[11]/n_ges_ch-P_ges_cap_t+x[15] #充电功率小于当前空闲
]
'''
第一行：CHP的功率情况
第二行：从气网买气的功率情况
第三行：EB的功率情况
第四行：P2G的功率情况
第五行：GB的功率情况
第六行：给电网卖电功率
第七行：从电网买电功率
第八行：EES充电功率
第九行：EES放电功率
第十行：TES充电功率
第十一行：TES放电功率
第十二行：GES充电功率
第十三行：GES放电功率
第十四行：EES储能状态
第十五行：TES储能状态
第十六行：GES储能状态
'''
#优化目标函数
def target_func(x):
  P_chp_t=x[0]
  P_ggrid_t=x[1]
  P_eb_t=x[2]
  P_p2g_t=x[3]
  P_gb_t=x[4]
  P_egrid_t_s=x[5]
  P_egrid_t_b=x[6]
  P_ees_t_cha=x[7]
  P_ees_t_dis=x[8]
  P_tes_t_cha=x[9]
  P_tes_t_dis=x[10]
  P_ges_t_cha=x[11]
  P_ges_t_dis=x[12]
  SOC_EES_T=x[13]
  SOC_TES_T=x[14]
  SOC_GES_T=x[15]

  # 计算电能NCI

  P_eload_t_net=P_eload_t
  a_ies_t_e=np.zeros(24)
  # for i in range(24):
  #   divided_sum=P_chp_t[i]+P_wt_t[i]+P_egrid_t_b[i]-P_eb_t[i]-P_p2g_t[i]
  #   if(divided_sum==0):
  #     a_ies_t_e[i]=0
  #   else: a_ies_t_e[i]=(a_chp_ec*P_chp_t[i]+a_grid_e*P_egrid_t_b[i])/divided_sum
  a_ies_t_e=(a_chp_ec*P_chp_t+a_grid_e*P_egrid_t_b)/(P_chp_t+P_wt_t+P_egrid_t_b-P_eb_t-P_p2g_t)
  a_pro_t_e=a_ies_t_e

  # 计算热能NCI
  a_ies_t_t=(n_t_t2e*a_chp_ec*P_chp_t+a_gb*P_gb_t)/(n_t_t2e*P_chp_t+n_t_e2t*P_eb_t+P_gb_t)
  a_pro_t_t=a_ies_t_t

  # 计算气能NCI
  a_ggrid_t_g=a_gas
  P_ggb_t=P_gb_t/(n_gb*LHV_ng)
  P_gchp_t=P_chp_t/(n_chp*LHV_ng)
  a_ies_t_g=(a_ggrid_t_g*P_ggrid_t)/(P_ggrid_t+n_p2g*P_p2g_t-P_ggb_t-P_gchp_t)
  a_pro_t_g=a_ies_t_g 

  # 计算电碳综合价格
  c_esp_t_ec=c_esp_t_e+k*a_pro_t_e
  # c_esp_t_ecc=c_esp_t_ee #此处为计算pro售电给ESP
  c_esp_t_tc=c_esp_t_t+k*a_pro_t_t
  c_esp_t_gc=c_esp_t_g+k*a_pro_t_g

  #目标函数参数
  I_egrid=np.sum((c_ebuy_t+k*a_ies_t_e)*P_egrid_t_s) #ESP出售给电网收益
  I_pro_e=np.sum(c_esp_t_ec*P_eload_t_net) #ESP出售电能给pro收益
  I_pro_t=np.sum(c_esp_t_tc*P_tload_t) #ESP出售热能收益
  I_pro_g=np.sum(c_esp_t_gc*P_gload_t)#ESP出售气能收益
  C_egrid=np.sum((c_esell_t+k*a_grid_e)*P_egrid_t_b) #购电支出
  C_ggrid=np.sum(c_gsell_t*P_ggrid_t) #天然气购买支出
  P_h_t=np.sum([P_chp_t,P_gb_t,P_eb_t,P_p2g_t,P_wt_t,P_pv_t],axis=1)
  C_op=np.sum(c_h*P_h_t) #设备维护支出
  C_en=k*c_gas*np.sum(P_gb_t+P_gchp_t) #环境支出


  return I_egrid+I_pro_e+I_pro_t+I_pro_g-C_egrid-C_ggrid-C_op-C_en


# 120000 120000
size_pop=60 #
dim1=16 #
dim2=24 #分为24个时间段
lb=np.zeros(24) #变量下限
ub=np.array([P_chp_max,P_ggrid_max,P_eb_max,P_p2g_max,P_gb_max,P_egrid_max,
             P_egrid_max,P_ees_max,P_ees_max,P_tes_max,P_tes_max,P_ges_max,P_ges_max,
             P_ees_cap,P_tes_cap,P_ges_cap]) #变量上限 16
max_iter=1500
prob_mut=0.3 #变异因子
CR=0.5 #交叉因子



de = MyDeFunc(obj_func=target_func,size_pop=size_pop,dim1=dim1,dim2=dim2,lb=lb,ub=ub,
              max_iter=max_iter,prob_mut=prob_mut,CR=CR,constraint_eq=constraint_eq,
              constraint_ueq=constraint_ueq)
de = de.run()

current_date = datetime.now().strftime('%Y-%m-%d_%H-%M-%S')
# 指定你想要创建的文件夹路径
folder_path = f'./fig/{current_date}'

# 检查文件夹是否已经存在，如果不存在则创建它
if not os.path.exists(folder_path):
    os.makedirs(folder_path)
    print(f"文件夹已创建：{folder_path}")
else:
    print(f"文件夹已存在：{folder_path}")

print("给电网卖电：",de.best_individual[5])
print("从电网买电：",de.best_individual[6])
print("EES充电：",de.best_individual[7])
print("EES放电：",de.best_individual[8])
print("TES充电：",de.best_individual[9])
print("TES放电：",de.best_individual[10])
# print("优化结果",de.population)
plt.rcParams['font.sans-serif']=['Microsoft YaHei']  # 指定默认字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号'-'显示为方块的问题




y=de.history_sum
plt.subplot(2,3,1)
plt.plot(y)
plt.grid(True)
# 添加标题和标签
plt.title('Profit')
plt.xlabel('迭代次数')
plt.ylabel('总利润')
# plt.show()
# file_name='Profit.png'
# full_path=os.path.join(folder_path,file_name)
# plt.savefig(full_path)
# plt.close

array = list(range(1, 25)) #横坐标




target_ees=de.best_individual[8]-de.best_individual[7]
target_tes=de.best_individual[10]-de.best_individual[9]
target_ges=de.best_individual[12]-de.best_individual[11]
target_egrid=de.best_individual[6]-de.best_individual[5]

plt.subplot(2,3,2)
plt.title('EES功率情况')
plt.xticks(array)
plt.xlabel('Time/h')
plt.ylabel('Power/kW')
plt.bar(array,target_ees,label='EES')

plt.subplot(2,3,3)
plt.title('TES功率情况')
plt.xticks(array)
plt.xlabel('Time/h')
plt.ylabel('Power/kW')
plt.bar(array,target_tes,label='TES')

plt.subplot(2,3,4)
plt.title('GES功率情况')
plt.xticks(array)
plt.xlabel('Time/h')
plt.ylabel('Power/kW')
plt.bar(array,target_ges,label='GES')

plt.subplot(2,3,5)
plt.title('Power Grid功率情况')
plt.xticks(array)
plt.xlabel('Time/h')
plt.ylabel('Power/kW')
plt.bar(array,target_egrid,label='Grid Power')

plt.subplot(2,3,6)
plt.title('EES电量')
plt.xticks(array)
plt.xlabel('Time/h')
plt.ylabel('Power/kW')
plt.bar(array,de.best_individual[13],label='EES CAP')
plt.show()
# plt.subplot(2,5,7)
# plt.title('TES电量')
# plt.xticks(array)
# plt.xlabel('Time/h')
# plt.ylabel('Power/kW')
# plt.bar(array,de.best_individual[14],label='TES CAP')

# 电能功率情况
# print("电能EES功率情况：",target_ees)
# print("电能Power Grid买电功率情况：",de.best_individual[10])
# print("电能Power Grid卖电功率情况：",de.best_individual[9])
# print("电能P2G功率情况：",de.best_individual[3])
# print("电能EB功率情况：",de.best_individual[2])
# print("电能CHP功率情况：",de.best_individual[0])
# bar2_=plt.bar(array,target_ees, label='EES')
# plt.subplot(2,5,8)
# bar1=plt.bar(array,target_ees, label='EES')
# # bar1_=plt.bar(array,target_ees, label='EES')
# bar2=plt.bar(array,target_egrid,color='red',bottom=target_ees, label='Power Grid')
# bar3=plt.bar(array,de.best_individual[3],color='green',bottom=target_egrid, label='P2G')
# bar3=plt.bar(array,de.best_individual[2],color='blue',bottom=de.best_individual[3], label='EB')
# bar3=plt.bar(array,de.best_individual[0],color='yellow',bottom=de.best_individual[2], label='CHP')
# plt.xticks(array)
# plt.xlabel('Time/h')
# plt.ylabel('Power/kW')
# plt.title('电能功率情况')
# plt.legend()
# plt.grid(True)
# # plt.show()
# # file_name='ElectricityPower.png'
# # full_path=os.path.join(folder_path,file_name)
# # plt.savefig(full_path)
# # plt.close()

# print("热能TES功率情况：",target_tes)
# print("热能GB功率情况：",de.best_individual[4])
# print("热能EB功率情况：",de.best_individual[2])
# print("热能CHP功率情况：",de.best_individual[0])
# # 热能功率情况
# plt.subplot(2,5,9)
# bar1=plt.bar(array,target_tes, label='TES')
# bar2=plt.bar(array,de.best_individual[4],color='red',bottom=target_tes, label='GB')
# bar3=plt.bar(array,de.best_individual[2],color='green',bottom=de.best_individual[4], label='EB')
# bar3=plt.bar(array,n_t_t2e*de.best_individual[0],color='blue',bottom=de.best_individual[2], label='CHP')
# plt.xticks(array)
# plt.xlabel('Time/h')
# plt.ylabel('Power/kW')
# plt.title('热能功率情况')
# plt.legend()
# plt.grid(True)
# # plt.show()
# # file_name='ThermalPower.png'
# # full_path=os.path.join(folder_path,file_name)
# # plt.savefig(full_path)
# # plt.close

# print("气能GES功率情况：",target_ges)
# print("气能Gas Grid功率情况：",de.best_individual[1])
# print("气能P2G功率情况：",de.best_individual[3])
# print("气能GB功率情况：",de.best_individual[4])
# print("气能CHP功率情况：",de.best_individual[0])
# # 气能功率情况
# plt.subplot(2,5,10)
# bar1=plt.bar(array,target_ges, label='GES')
# bar2=plt.bar(array,de.best_individual[1],color='red',bottom=target_ges, label='Gas Grid')
# bar3=plt.bar(array,de.best_individual[3],color='green',bottom=de.best_individual[1], label='P2G')   #需修改
# bar3=plt.bar(array,de.best_individual[4],color='blue',bottom=de.best_individual[3], label='GB')
# bar3=plt.bar(array,de.best_individual[0],color='yellow',bottom=de.best_individual[4], label='CHP')
# plt.xticks(array)
# plt.xlabel('Time/h')
# plt.ylabel('Power/kW')
# plt.title('气能功率情况')
# plt.legend()
# plt.grid(True)
# plt.show()
# file_name='ESP.png'
# full_path=os.path.join(folder_path,file_name)
# plt.savefig(full_path,dpi=600)

target_func(de.best_individual)





























'''
def obj_func(x):
    return np.sum(x)

constraint_eq=[
    # 等式约束，默认 = 0
    lambda x:8 - np.sum(x[0]),
    lambda x:12 - np.sum(x[1])
]

constraint_ueq=[
    # 不等式约束，默认 <= 0
    lambda x: x[0][0] + x[0][1] -5,

]


[
  [[1,2,3],[1,2,3]],
  [[1,2,3],[1,2,3]]
]



# 假设每个变量数组有两个元素
iter=500 #迭代次数
size_pop = 60  # 种群大小
dim1 = 2  # 有几个变量
dim2 = 3 # 每个变量状态分为24个小时
lb = np.zeros(2)  # 每个元素的下界
ub = np.array([3,4])  # 每个元素的上界
nprocess=0 #是否多线程
pro_mut=0.5 # 变异因子
CR=0.5





# 初始化种群
# population = np.random.uniform(low=np.array(lb), high=np.array(ub), size=(size_pop, dim))
de = MyDeFunc(obj_func=obj_func,size_pop=size_pop,dim1=dim1,dim2=dim2,
              lb=lb,ub=ub,max_iter=iter,prob_mut=pro_mut,CR=CR,
              constraint_eq=constraint_eq,constraint_ueq=constraint_ueq)
print('population_raw:\n:',de.population)
de = de.run()
print('population_fitness:\n:',de.population)
x=np.random.randint(0,iter+1)
y=de.history_sum
plt.plot(y)
plt.grid(True)
# 添加标题和标签
plt.title('Line Graph Example')
plt.xlabel('x')
plt.ylabel('y')

plt.show()
'''